%%
%Code for creating the main part and marginals of Figure 4
load('example_FOV_data.mat')

%%%%%%%%%%%%%%%%%
%Next pick a day and plot out cells color-coded by their responses
days = [2 6 14];
mouse = 1;
cmap = hot;
cmap = cmap([1:3:64],:); cmap = flipud(cmap); widths = [1:3:64]./12;
nameDays = {'Day -1'; 'Day 1'; 'Day 18'};
dBs = [20:10:80];

x = linspace(0,512);
y = linspace(0,512);
[XX,YY] = meshgrid(x,y);
pred = [XX(:),YY(:)];
p = predict(mdl,pred);
f = @(x) -(x*mdl.Beta(1) + mdl.Bias)/mdl.Beta(2);
y = f(x);



figure
hold on
for k = 1:3


for intensity = [1:7]
    
day = days(k);
subplot(3,7,((k-1)*7+intensity))

for i = 1:length(allFuncs{mouse,day}) %Loop through cells
    
    %Get boundaries
    bounds = allData{day}.cellInfo.dat(i).boundaries{:};
    
    %Use response to get color (for certain intensity)
    if isResp{mouse,day}(i) %If responsive to this freq
        
        resp = allFuncs{mouse,day}(i,:); resp = resp(intensity); %Get response to this intensity
        color_idx = round(resp);
        if color_idx < 2; color_idx = 1;  end %Account for small responses
        if color_idx > length(cmap); color_idx = length(cmap); end %Account for large responses
        color = cmap(color_idx,:);
        if color_idx == 1; color = [0 0 1]; end
    else
        %For non-responsive cells just use the darkest color
%         color_idx = 1;
%         color = cmap(color_idx,:);
        color = [0 0 1]; color_idx = 1;
        
    end
    
    if color_idx ~= 1
        for j = 1:length(bounds) %For each entry

            %Account for registration offset
            bounds(j,1) = bounds(j,1) + allData{day}.cellInfo.yrange(1) - 1;
            bounds(j,2) = bounds(j,2) + allData{day}.cellInfo.xrange(2) - 1;

        end
   
   
    plot(bounds(:,1),bounds(:,2),'linewidth',widths(color_idx),'color',color);
    hold on
    end
    
end
set(gca,'xtick',[])
set(gca,'ytick',[])
box off
plot(y,x,'--','Color',[0.3 0.3 0.3],'LineWidth',2.5)
xlim([0 512])
ylim([0 512])

if intensity == 1
    ylabel(nameDays{k});  
end
if k == 1 
    title(strcat(num2str(dBs(intensity)),' dB SPL')); 
end

end

end

%%%%%%%%%%%%%%%%%%%%%%%
%Now create the marginals. First start with the tuning curves

freq = 8000;
days = [2 6 14];
cmap = hot;
cmap = cmap([1:3:64],:); cmap = flipud(cmap); widths = [1:3:64]./50;
mouse = 1;
intensities = 20:10:80;
day_names = {'Day -1' 'Day 1' 'Day 18'};

for k = 1:3 %Loop through days
    
day = days(k);
figure
hold on

for i = 1:length(allFuncs{mouse,day}) %Loop through cells
    
    %Use response to get color (for certain intensity)
    if isResp{mouse,day}(i) %If responsive to this freq
        
        resp = allFuncs{mouse,day}(i,:); resp = max(resp); %Get response to this intensity
        color_idx = round(resp);
        if color_idx < 2; color_idx = 1;  end %Account for small responses
        if color_idx > length(cmap); color_idx = length(cmap); end %Account for large responses
        color = cmap(color_idx,:);
        if color_idx == 1; color = [0 0 1]; end
    else
        %For non-responsive cells just use the darkest color
        color = [0 0 1]; color_idx = 1;
        
    end
    
    func = allFuncs{mouse,day}(i,:);
    
    if color_idx ~= 1
        
        if any(ismember(find(func == min(func)),[1:(round(length(func)/2)+1)])) && any(ismember(find(func == max(func)),[(round(length(func)/2)+1):length(func)])) %Check if it's monotonic-ish
      
        plot(intensities, func,'linewidth',widths(color_idx),'color',color);
        hold on
        
        end
    end
    
end

ylabel('Z-Scored Spikes')
xlabel('Intensity')
box off
ylim([-5 100])
set(gcf, 'Position', [100 100 500 600])
title(day_names{k})

end


%%%%%%%%%%%%%%
%For making the distance marginals
days = [2 6 14];
intensity = 6;
mouse = 1;
nGroups = 5;
colors = [0.8 0.8 0.8; 0.5 0.5 0.5; 0.2 0.2 0.2]; 

figure
for intensity = 3:7

subplot(1,5,(intensity-2))
for k = 1:3 %Loop through days
    
clearvars dists resps gDists gResps
day = days(k);
count = 1;

for i = 1:length(allFuncs{mouse,day}) %Loop through cells
    if isResp{mouse,day}(i) %If responsive to this freq
        nums = allFuncs{mouse,day}(i,:); resps(count,1) = nums(intensity);
        dists(count,1) = distances{mouse,day}(i,:);
        count = count + 1;
    end
end

%Now sort the data
[dists,inds] = sort(dists); resps = resps(inds);
    
%Next go through the neurons and get average distance and response measures
num1 = ceil(length(resps)/nGroups); num2 = mod(length(resps), nGroups);
groups = ones(1,nGroups).*num1; 
if num2 ~= 0; groups(1:(nGroups-num2)) = groups(1:(nGroups-num2)) - 1; end

start = 0; stop = 0;
for i = 1:nGroups
    start = stop + 1; stop = stop + groups(i);
    gResps(i) = mean(resps(start:stop)); gDists(i) = mean(dists(start:stop)); 
    gError(i) = std(resps(start:stop))/sqrt((stop - start));
end


plot(smooth(dists,50),smooth(resps,50),'Color',colors(k,:),'LineWidth',2)
hold on
end
legend({'Day -1' 'Day 1' 'Day 18'})
xlabel('Distance from SVM Line (microns)')
ylabel('Average Response')
xlim([-300 500])
ylim([0 35])
box off

end